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Abstract 

We used A^-body dynamical simulations to analyze the conditions for 
the gravitational stability of a three-dimensional stellar disk in the grav- 
itational field of two rigid spherical components - - a bulge and a halo, 
whose central concentrations and relative masses vary over wide ranges. 
The number of point masses N in the simulations varied from 40 to 500 
thousands and the evolution of the simulated models is followed over 10- 
20 rotation periods of the outer edge of the disk. The initially unstable 
disks are heated and, as a rule, reach a quasi-stationary equilibrium with 
a steady-state radial-velocity dispersion c r over five to eight periods of 
rotation. The radial behavior of the Toomre stability parameter Qt{t) 
for the final state of the disk is estimated. Numerical models are used 
to analyze the dependence of the gravitational stability of the disk on 
the relative masses of the spherical components, disk thickness, degree 
of differential rotation, and initial state of the disk. Formal application 
of existing analytical local criteria for marginal stability of the disk can 
lead to errors in c r of more than a factor of 1.5. It is suggested that the 
approximate constancy of Qt — 1.2 1.5 for r ~ (1-^2) x L (where 
L is the radial scale of a disk surface density), valid for a wide range 
of models, can be used to estimate upper limits for the disc mass and 
density based on the observed distributions of the rotational velocity of 
the gaseous component and of the stellar velocity dispersion. 



2 



1 INTRODUCTION 

Galactic disks, which consist mostly of old stars, can be considered as collisionless systems in 
quasi-stationary equilibrium with very slow evolution. This state is characterized by certain 
radial dependences of the stellar velocity dispersions (c r , c^,, c z ) that ensure stability of the disk. 
Knowledge of the stability conditions makes it possible to develop self-consistent models for the 
disks of real galaxies for which both the rotational velocities and stellar velocity dispersions have 
been measured at various galactocentric distances. 

The problem of determining the minimum stellar velocity dispersion sufficient to ensure sta- 
bility of the disk against arbitrary perturbations is especially important if, as many authors have 
suggested (see, e.g., [1-7]), real galaxies may be in a state of threshold (marginal) stability. This 
approach enables the local density or integrated mass of the disk to be estimated from the ob- 
served velocity dispersion. In the general case, the old stellar population of a galactic disk can 
have an excess velocity dispersion in the presence of other factors that heat the disk, which are 
not directly related to the gravitational instability. However, even in this case, the conditions for 
marginal stability provide valuable information by yielding an upper limit for the mass of the disk 
that enables it to be stable. 

Fridman and Polyachenko [8, 9] carried out a detailed theoretical analysis of the stability of 
thin rotating disks against different kinds of perturbations (including bending modes). Together 
with certain advantages over numerical simulations (the mathematical rigorousness of the solu- 
tions in the framework of the problem formulated), the analytical approach to the dynamics of 
perturbations in a gravitating disk and the conditions for 

stability has the drawback that it can be implemented only for very simple models and can 
yield only coarse estimates for the parameters of the disk component when applied to real objects. 
Let us consider these simplifications in more detail. 

First and foremost, the simple analytical models of collisionless disks that are used in stability 
analyses usually assume that the disk thickness is small; they actually consider an infinitely thin 
layer. Despite the smallness of the ratios h z /r and h z /L over most of the stellar disk (here, r 
is the radial coordinate, h z is the vertical scale height, and L is the radial scale length of the 
disk), this condition may not be sufficient to justify neglect of vertical motions (see Gor'kavyi and 
Fridman [10, Appendix II 1 ] for a detailed discussion of this issue; in the general case, the dynamics 
equations for astrophysical disks cannot be adequately treated in a two-dimensional formulation). 

Analytical studies of the dynamics of small perturbations in a stellar disk usually assume that 
the disk parameters have small radial gradients (see, e.g., [1, 11, 12]). Attempts to allow for 
gradients have thus far been made only in terms of WKB approximations (see, e.g., [2, 13-15]). 
This means that the characteristic perturbation wavelength is short compared to the local scale 
lengths for variations of the circular velocity V(r), stellar radial-velocity dispersion c r (r), and 
surface 

density cr(r). In many cases, these conditions are barely satisfied or even not satisfied. Taking 
into account differential rotation may pose the most difficult problem. Allowance for weakly dif- 
ferential rotation is possible for small nonradial perturbations [13, 16]. However, nonaxisymmetric 

1 Written in coauthorship with O.V. Khoruzhii. 
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(in the limiting case, spoke-like) perturbations are more unstable. Their stabilization under the 
same conditions requires significantly stronger disk heating, i.e., a higher stellar velocity disper- 
sion [17], and these very perturbations apparently increase the stellar velocity dispersion at the 
nonlinear stage in the case of initially cool systems. 

Analytical studies of the gravitational stability of disks usually apply the epicyclic approxi- 
mation c r <C V, which is valid at the peripheries of most galaxies, where c r /V ~ 0.1 4- 0.3, but 
breaks down near their centers, where c r >V. 

Another important limitation of analytical approaches is that the criteria derived are local, 
whereas a number of studies suggest that the disk stability conditions have a global nature [18, 
19]. This means that the equilibrium parameters, e.g., at the center, while leaving their values at 
the disk periphery unchanged, may affect the conditions for gravitational stability throughout the 
disk. In contrast to a local approach based on the analysis of dispersion equations, global analyses 
aim to determine the eigenfrequency for the entire disk by solving a boundary-value problem, 
which determines the influence of the conditions in one part of the system on the dispersion 
properties of the perturbations throughout the disk. Therefore, in a rigorous approach, the disk 
must be considered as a whole. However, global analyses have been performed only for certain 
specific power-law distributions [18-20]. For example, Bertin et al. [21] considered global modes 
in galactic disks as possible mechanisms for maintaining long-lived spiral density waves. 

Numerical simulations of collisionless systems are more flexible in terms of the choice of model. 
They make it possible to go beyond simple two-dimensional models and directly follow the de- 
velopment of perturbations in a disk that is initially in equilibrium. However, this approach has 
drawbacks of its own. The most serious problems of iV-body simulations include (1) certain math- 
ematical simplifications that are inevitable when a disk is modeled as a system of N gravitating 
bodies, where N is incomparably smaller than the number of stars in real galaxies, and (2) the 
dependence of the final state of the system (after it reaches quasi-equilibrium) on the initial pa- 
rameters, which are poorly known for real galaxies. When comparing simulation results with real 
galaxies, it can also be difficult to allow for the dissipative galactic medium (gas), in which the 
sound speed is much lower than the stellar velocity dispersion. We consider these problems in 
more detail below. 

The principal goal of this study is to determine for galaxies with various mass distributions 
the minimum local disk velocity dispersions that enable their three-dimensional disks, which are 
initially in a weakly unstable equilibrium, to reach a quasi-stationary state. Sections 2 and 3 
give a concise review of previous results of analytical and numerical approaches to estimating the 
stellar velocity dispersions required to ensure the stability of collisionless disks. Section 4 describes 
the principles underlying the construction of the dynamical models and the determination of the 
stability threshold used in this paper. Section 5 considers various three-dimensional models of 
galaxies, and the last section presents and discusses the main results. 
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2 ANALYTICAL STABILITY CRITERIA 

Several criteria for gravitational instability derived analytically using various models have been 
published. Let us review those we consider to be most important. 

(a) The Toomre Criterion 

In order for an infinitely thin, uniform, rigidly rotating stellar disk to be gravitationally stable 
against axisymmetric perturbations, it must satisfy the following condition derived by Toomre [1] : 

3.36G(T c r 
c r >c T = , Qt = — > 1, (1) 

* c T 

where 6e = 2Q^Jl + rdQ/2Qdr is the epicyclic frequency and a is the surface density. Condition (1) 
assumes that the epicyclic approximation is valid, i.e., that the difference between the velocity 
V(r) and the circular velocity V c (r) can be neglected. Although this relation was derived in a 
local analysis, the study of Evans and Read [19] of the eigenmodes for self-similar disks in a 
corresponding approximation overall supports its validity. Miller [22] compared the theoretical 
increments derived for a Toomre model with the results of simulations of axisymmetric modes in 
which all other perturbations were artificially suppressed. Experimental gradients were shown to 
be consistent with the theoretical results. 

(b) Allowance for Finite Disk Thickness 

Finite thickness is a stabilizing factor for gravitational instability in the plane of the disk [9, 
12]. The generalization of the Toomre stability criterion (1) for the case of a disk of finite thickness 
obtained by Morozov [13, 14, 16] has the form 

Qt ] = < 1, (2) 

Vt 1 + 0.974Aa;/c r ' 1 } 

where A is the half-thickness of the isothermal self-gravitating disk. However, this condition 
proved to be far from sufficient to ensure the stability of real systems. This became clear from 
the first numerical simulations performed in the 1970s, which showed that c r ~ (1.5 -j- 5) x ct at 
the periphery of a stationary, collisionless disk [23-30]. 

(c) Simplified Allowance for Nonaxisymmetric Perturbations 

The stronger instability of spiral waves compared to axisymmetric modes is one of the factors 
that makes the Toomre criterion inadequate. Polyachenko and Shukhman [31], Kalnajs [32], and 
Polyachenko and Fridman [9] were the first to show that nonaxisymmetric modes are the dominant 
instabilities in a gravitationally unstable disk. As pointed out by Polyachenko and Fridman [9], the 
azimuthal- velocity dispersion c v is smaller than c r (except for the innermost regions). Therefore, 
the relation 

Clp = Cr 2ir (3) 

implies that, in order for an azimuthally cooler disk to be stabilized, it must be heated more 
strongly, so that, in view of (3), the condition (1) can be rewritten in the form 

Q { ? > S, where S = — . (4) 
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The parameter S characterizes the degree of differential rotation of the disk. In most cases, we 
can assume that 1 < S < 2, based on the observed shapes of galactic rotation curves. This form 
of the stability condition is discussed in [13-16]. 

Thus, the azimuthal velocity dispersion determines the elasticity of the medium against strongly 
nonaxisymmetric perturbations, so that the suppression of gravitational instability requires stronger 
heating of the disk by a factor of 20,/ se.. 

Criterion (4) can be considered to be the Toomre criterion with a simplified allowance for 
nonaxisymmetric perturbations, 
(d) The Morozov Criterion 

Morozov [13], Morozov and Khoperskov [14], and Morozov [16] analyzed the dynamics of 
weakly nonradial perturbations in a nonuniform disk in a WKB approximation. The resulting 
stability criterion takes into account many factors (radial nonuniformity of the surface density a 
and radial- velocity dispersion c r , the disk thickness, and differential rotation): 



Q ( T M) 
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(5) 



where S = 2Q/ae and D = (1 + 0.974a;A/( 1 S 2 c T ))- 1 . 

In the general case, determining c^ M ) from (5) amounts to integrating a reduced differential 
equation [33]. The main drawback of this criterion is that it was derived in the context of the 
dynamics of tightly wound spiral waves (m/r <C k r , where k r is the radial wavenumber) and then 
formally applied to spokelike perturbations. 

If the scale length L c = (dlnc^ j 'dr)~ l is fixed, the differential equation (5) reduces to a 
simple algebraic relation [13, 14, 16]. For typical rotation-curve shapes, we have D ~ 0.6 -r- 0.8. 
Introducing the factor D also makes it possible to formally allow for the disk thickness in criteria (3) 
and (6). 

(e) The Polyachenko— Polyachenko— Strel'nikov Criterion 

Unlike the criteria discussed above, the analysis of Polyachenko et al. [17] focused on extremely 
nonaxisymmetric perturbations in a thin disk. Under the assumptions made, the rotation curve 
is the sole factor determining the stability boundary, since this boundary depends only on the 
parameter n = —rdQ/(Qdr). Figure 1 in [17] shows the dependence of the dimensionless velocity 
dispersion at the stability boundary on the parameter a = 2/n. We will use the approximating 
function 

(6) 



Qt ^ ~ 



u 



'1.1 



c T V exp(a — 1/4) — 1 ' 

where is the minimum Toomre parameter for a stable disk. 

This approximation has sufficient accuracy for our needs when 1.2 < a 2 . The criterion of 

Polyachenko et al. [17] depends on a single parameter, since the form of the rotation curve fully 

(p) 

determines the stability boundary. In particular, it follows that Q T ~ 3 in the region with 
constant rotational velocity (n = 1). 

(f) Allowance for the Gaseous Subsystem in a Stability Analysis for Stellar Disks 

The presence of a cooler component also contributes to destabilization of the stellar disk. A 
number of authors [34-38] have analyzed this problem in detail as applied to radial perturbations. 
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Ortega et al. [38] considered a more general problem: they analyzed how the inhomogeneous 
composition of a thin disk influences the stability of the disk against small radial perturbations 
when the disk consists of particles having a mass spectrum such that more massive particles have 
lower velocity dispersions. 

The problem is simplified if the mass of the "cool" component is relatively small. Assuming 
the gas surface density <7gas to be usually much lower than the surface density of the stellar disk 
a star an d c s <S c r (where c s is the adiabatic sound speed in the gas), onecan write for the velocity 
dispersion of a stellar disk containing gas at the stability boundary [37] 

ct o-g as + a star 1 + (c s /c T ) 2 

For example, for the parameters of the solar neighborhood in the Galaxy, it follows from (7) that 
the gaseous subsystem, which contributes about 20% of the disk surface density, increases the 
minimum stellar-velocity dispersion sufficient for stabilizing radial perturbations by < 10%. 

Finally, we note that attempts have also been made to determine the parameters of the disk 
subsystem using other approaches based on the possible existence of various structures in the 
disk (spiral density waves or a bar) rather than on local stability conditions, assuming certain 
mechanisms that form and sustain these structures. However, the analysis of these approaches is 
beyond the scope of this paper. 



3 ANALYSIS OF GRAVITATIONAL STABILITY IN NU- 
MERICAL SIMULATIONS 

Numerical simulations describing the dynamical evolution of disks can be used to analyze disk 
instabilities for given initial conditions allowing for nonuniform distributions of the mass and an- 
gular velocity. Such simulations naturally take into account the formation of appreciably nonlinear 
and nonaxisymmetric structures such as bars or transient spirals. The main problem with this 
approach is that the results depend on the initial conditions, since the evolution of a disk starting 
from a strongly unstable state can proceed very differently from the evolution of a disk starting in 
a subcritical state. In addition, the disk can suffer from slow secular instabilities that are difficult 
to take into account in numerical simulations. 

N-bodj evolutionary models are usually aimed at studying the development of instability or 
- which is of the most interest for us — at analyzing the states of a system after many disk 
rotations. Of the large number of published studies of this type, we will mention those that are, 
in our view, most important in the context of establishing the conditions for stability of a disk. 

Carlberg and Sellwood [39, 40] analyzed the influence of small, nonstationary perturbations 
of the potential on the evolution of the velocity-distribution function. In particular, they analyzed 
scattering by nonstationary spiral waves; the resulting growth of the stellar-velocity dispersion 
agrees well with the results of numerical simulations. 

The critical velocity dispersion for stellar disks (the parameter Qt) has been computed many 
times based on the results of iV-body dynamical simulations (see, e.g., [7, 24, 28, 33, 40-46]). None 
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of the simulations yielded a stable disk with < 1. As a rule, QT{ r ) increases with distance from 
the center of a galaxy. This pattern appears both in simulations of three-dimensional disks [5, 
46-48] and in some theoretical analyses [49] . 

The work of Athanassoula and Sellwood [34] can be considered to be classic. They concluded 
that a two-dimensional disk is always stable if its Toomre parameter exceeds Qt >2.2 + 2.4. 
However, this conclusion did not take into account vertical motions. Since three-dimensional disks 
are gravitationally more stable, two-dimensional models underestimate the mass of a marginally 
stable disk. In addition, Athanassoula and Sellwood [34] used a radially averaged Qt and simulated 
a specific density distribution based on a Toomre-Kuz'min model. It is therefore not clear whether 
their results are applicable for real three-dimensional disks. 



4 NUMERICAL MODELS: SPECIFYING INITIAL CON- 
DITIONS 

The dynamical models used in this paper are based on numerical integration of the equations of 
motion for N gravitationally interacting particles taking into account the external field produced 
by the steady-state mass distribution in the bulge and halo. 

For the halo, we adopted the commonly used spatial distribution of the volume density in the 
form 

gfe(r)= ( i+eW ' (8) 

where £ = Vr 2 + z 2 is the radial coordinate. Choosing k = 1 yields a flat rotation curve in the 
halo-dominated region. The relative central density of the spherical component increases with k, 
imitating the influence of the bulge. Because the bulge enters our model as a separate spherical 
component, we restrict our analysis to a halo model with k — 1. 

The following set of equations describes the dynamics of the N gravitating bodies: 

^ = J2h + F s (i = l,...,JV), (9) 

— # 

Here, the radius vector fi(t) determines the position of the ith particle, is the force of interaction 
between the ith and jth particles, and the force F s = Ff, + is due to the bulge/halo spheroidal 
subsystem. The halo mass distribution (8) with a central density of Qho = Mh/{4:ira 3 [R/a — 
tan(i?/a)]} yields for the force 

which is determined by the spatial scale lengths a and mass Mh inside the sphere £ = |r| < R. 
We adopt a King model for the density distribution in the spherical bulge: 

„ f Qbo/ [1 + {i/b) 2 f 2 , £<(r b ) max ( , 
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where 

M b = Anb 3 gJ In \(r b ) max /b + yl + ((r 6 ) max /6) 2 l - foW& = | (12) 

1 V 1 + (( r i)max/6) 2J 

is the mass of the bulge. We have for the gravitational force inside £ < (r 6 ) ma x 



It is obvious that F b = —GM b r/£ 3 in the domain r > (r b ) ma x- 

The dynamical model must adequately describe the Newtonian interactions between stars and 
ensure that the system is collisionless. This is achieved by modifying the gravitational force at 
small distances by introducing a potential-cutoff radius r c for any pair of interacting particles % 
and j. The optimum choice of cutoff radius and number of particles has been widely discussed in 
the literature (see, e.g., [50-53] and references therein). 

Here, we use a Plummer model for the potential: 

% = ~G ! ; \ , (14) 

where is the distance between particles and r c is the cutoff radius. For a fixed number of 
particles N, it is always possible to choose a cutoff radius r c that ensures that the model is 
collisionless. However, the number of particles must be sufficiently large to minimize the error 
introduced by the modification of the particle-interaction potential at small distances. 

We characterize the disk surface density by the scale length L, which determines the exponential 
law <r(r) = cr exp(— r/L). We assumed that the disk surface density is zero in the region r > 5L 
at the beginning of the simulations. We used a system of units in which G — 1, R — AL — 1, and 
the mass of the disk is = 1. We normalized the mass of the halo Mh inside the radius £ < AL 
to the mass of the disk, \i = M h jM^. In this system of units, one period of rotation of the outer 
edge of the disk lies in the range t ~ 3 4- 4 for /j, = 1 4. 

Let us now describe the procedure for specifying the initial density distribution along z axis 
and the residual velocities of the equilibrium disk. 

The vertical structure of the disk is determined by the equations [54] 

1 d ( <9$\ <9 2 $ A . 2 dg 9$ 

-rlfr {'IF ) + fl? = 4 ^ + ' C ^ = ~Tz ^ (15) 

where g and ^ s are the spatial density in the disk and spheroidal subsystems, respectively, and c z is 
the vertical-velocity dispersion, which is assumed to remain constant with z at t — 0. Eliminating 
the potential $ from (15) and introducing the circular velocity 



V c (r) = 



\ 




r Ur , (16) 



z=0 



we can transform the set of partial differential equations to an approximate equation for the disk 
volume density g(z) in the form of an ordinary differential equation [54]: 
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Together with the conditions g(z = 0) = go, dg(0)/dz = 0, and / g(z; r)dz = cr(r), this equation 

— oo 

determines the vertical structure of the disk at the radius r for a given surface-density distribution 
a. The E term can produce a large error in the estimated density at the disk center, and, in 
practice, it is assumed that E(r — > 0) — > 0. To determine go and g(z) for specified g s (z, r), c z (z, r), 

oo 

and cr(r), we construct the function F(g Q ) — 2 / g(z)dz — a. We solve the equation F(qo) = 

o 

iteratively jointly with numerical integration (17). After determining the density distribution in 
z, the particles are arranged along the vertical axis on a grid Zk = kAz (k = —K, . . . , K). The 

fcth cell contains particles in proportion to Ofc/cr, where = J g(z)dz. 

Zk-l 

Equation (17) is approximate, since it was derived neglecting the dependence of the potential 
on the vertical coordinate in the first term in (15). Strictly speaking, a disk constructed in this 
way is not in equilibrium. However, we are interested in initial, unstable states that evolve to new 
stationary states of the disk. Therefore, the lack of an exact equilibrium in the vertical direction 
plays the role of a small additional initial perturbation. 

The initial velocity distribution is a Schwarzschild function and has the form of an anisotropic 
Maxwellian distribution: 

„. ( u 2 (v - rfi) 2 w 2 } 

t (u,v,w) = A exp < > , 

I 2c 2 2c% 2c 2 z j ' 

where {u,v,w} are the velocity components of the particles in cylindrical coordinates. 

To obtain the model with the minimum velocity dispersion in the final, stable disk state, 
we chose a subcritical disk state whose dynamical evolution ensured both the preservation of 
an exponential surface-density profile and stability against initial perturbations in the plane of 
the disk and against bending perturbations (responsible for the increase in the vertical velocity 
dispersion) at the end of the computations. 

In practice, the initial radial-velocity dispersions c r in models with low-mass bulges corre- 
sponded to Toomre parameters Q T — 0.8 + 1.1 and Q T — 1-2 4- 2.2 in the central region (r <2L) 
and at edge of the disk, respectively. In the case of massive bulges, the models started from higher 
central values of Qt- The initial vertical-velocity dispersion was set c z proportional to c r . We 
considered two types of models with different initial disk thicknesses: "thin" disks unstable in the 
z direction 2 , which had central (c z /c r )o values (c z /c r )o = (0.4 0.5) at t — 0, and "thick" disks 
with (c z /c r ) = (0.6 -i- 0.8), which are close to the stability limit against bending perturbations. 
We assumed that these ratios varied slowly in radius in accordance with an exponential law with 
a radial scale length appreciably exceeding L. This choice of velocity dispersion ensured weak 
instability of the disk at all r. With the exception of specially stipulated cases, we describe below 
numerical models of "thick" disks. 

In none of the models did the velocity dispersion c r remain constant: the disk "heated up" in 
the course of its evolution. 

2 The bending instabilities of the oscillation mode lead to the heating of a thin disk in the vertical direction and, 
as a result, to its thickening. As our numerical models show, the axisymmetric bending oscillation mode plays an 
important role. 
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We found the mean tangential velocity of the model point masses by solving the Jeans equation 
assuming the absence of systematic radial motions, axial symmetry, and symmetry about the plane 
z = 0: 

1 ; c r \ cy gcj dr c 2 r dz J ' K ' 

Here, < . . . > denotes averaging of the velocities and the last term in (18) is due to the chaotic 
components of the radial velocity u and vertical velocity w. When specifying the initial state of 
the model, we assumed that < u >= and < w >= and assigned the rotational velocity of the 
disk in accordance with (18). Thus, initially we have a balance of the radial and vertical forces, 
and hence the disk begins to evolve from a nearly equilibrium state. 

In all the computations, we specified the initial distribution of the velocity dispersion c v at 
t = in accordance with (3). We verified that the condition Q c = — — - = 1 was satisfied as the 

disk evolved. In models with fairly massive halos (// > 2), the mean deviations of Q c from unit at a 
given r did not exceed 3% over several dozen rotations of the disk edge. This error is partially due 
to the numerical differentiation used in a process of computing the epicyclic frequency as. In the 
models with massive halos, the domain in which c r jse > r is small (for /i — 3, we have r <0.03). 
In the case of // <1, this domain expands to r <0.15. In addition, the vertical disk scale length 
increases in such models, and these factors result in stronger deviations from the equality (3). The 
amplitude of fluctuations Q c {t) decreases as the number of particles increases and does not exceed 
2% in models with jj, >2 and iV = 2 x 10 5 (except for the innermost region of the disk, r <0.5L), 
which reflects the collisionless nature of the models constructed. 

Using the approach described above, we performed more than 40 numerical simulations of the 
dynamical evolution of a disk toward a steady state in dependence of the initial velocity dispersions 
c r (r) and c z (r). We considered a wide range of parameters for the bulge (M b /M^ = 3, 
b/L = 0.02 + 0.8) and halo (M h /M d = 3.5, a/L = 1 4). The number of particles was 
N = (40 + 500) x 10 3 in the TREEcode computations. 

We verified the stability of the solutions against the choice of computational method by com- 
paring results for several models obtained using two very different methods to compute the gravi- 
tational force: the TREEcode method and direct "particle to particle" (PP) integration, in which 
each particle interacts with each of the other particles, for iV = (20 -r 80) x 10 3 . As an example, 
Figure 1 shows the time dependence of the radial-velocity dispersions for one of the models com- 
puted using these two methods, with similar initial conditions. 3 A comparison of the two results 
reveals no significant differences between the final disk states. 

3 The velocity-dispersion components in the plane of the disk c r and c v increase rapidly during the initial stage 
(t < 10) of the heating of an initially cool disk with initial Toomre parameter Qt{t < 2£) ~ 0.85. The relaxation 
of the disk in the vertical direction is appreciably slower, so that small local decreases of c r are possible when the 
disk is heated in z or there are radial motions (see curve 1 in Fig. 7 after t > 7). 
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5 DETERMINING THE THRESHOLD Q T VALUES 

5.1 Disk Heating Mechanism 

An initially axisymmetric, equilibrium disk is heated, increasing its velocity dispersion with time. 
The question of greatest importance for dynamical simulations is the heating mechanism. To 
correctly describe the processes in stellar disks, it is important, in particular, to ensure that the 
heating is not due to the collisional relaxation of the particles, whose number iV is many orders 
of magnitude smaller than the number of stars in real 

systems. This is achieved via appropriate choices of the number of particles and the cutoff 
radius. Our criterion for the absence of significant collisional-relaxation effects is that the character 
of the system's evolution be preserved when the computations are performed for an increased 
number of particles, as we verified for many models 4 . 

Below, we summarize the most important features of the disk heating shown by the numerical 
simulations. 

(1) The disk-heating time is much greater than the mean rotation period of the particles 
(Fig. 2a). In the initial stage (t <1), c r remains virtually constant, as long as the disk remains 
axisymmetric. In the case of a low-mass or nonexistent halo, the evolution of the disk is determined 
by the bar mode, and the disk is heated due to the formation of an nonaxisymmetric bar and 
associated two-armed spiral. Models with sufficiently massive haloes do not show any enhancement 
of the bar mode, however they develop a complex transient system of small-scale spiral waves 
(Fig. 2b). The decrease of the amplitudes of these waves with time is accompanied by the transfer 
of rotational kinetic energy to the chaotic component of the velocity, resulting in heating of the 
disk. 

(2) The heating of an initially cool disk (0.5 < Qt ^ 1) begins in its central region (Figs. 2a, 2b). 
The heating at the periphery proceeds much more slowly. The velocity dispersion at the periphery 
usually begins to rise when the center has already reached a quasi-stationary state (Fig. 2a). At 
the same time, the processes at the center and at the periphery are interrelated: fast growth 
of instability in the central region can speed up the heating of the outer part of the disk, while 
stability of the central region can slow down this process. 

3) If the system does not develop a bar, the amplitude of the perturbations begins to decrease 
with increasing velocity dispersion. In turn, the increase of the radial-velocity dispersion c r slows 
down parallel with decreasing wave amplitude. The heating virtually ceases after the decay of the 
transient spiral waves (see Figs. 2a, 2b). The integrated amplitude of the Fourier harmonics 

/ 1 N 

A(m;t) = l^2A 2 (m,p;t) , A(m,p,t) = — J2 ex P{* [ m ^j(0 +P m ( r j(0)]} 

4 The minimum number of particles N cr ^ for which this criterion is satisfied to sufficient accuracy (i.e., to within 
random fluctuations of the estimates of the final parameters) depends, in particular, on the form of the rotation 
curve. The number N CT ^ is larger for highly concentrated nuclei (bulges), due to the higher rotational angular- 
velocity gradient in the central region of the disk. In the limiting case of a very short, rigidly rotating region, N CT ^ 
can reach 3 x 10 6 [20]. Our simulations showed that, if the rotation curve grows monotonically to r ~ 2L, then 
-^crit can k e se t equal to ~ 4 x 10 4 . 



12 



decreases with time for all mode numbers m, but most slowly for m = 2 (Fig. 2c). The density 
distribution in the disk becomes almost axisymmetric (if the mass of the spherical components is 
large enough to prevent the development of a bar), and, on the whole, the velocity dispersion c r 
maintains its level over several dozen rotations, provided that relaxation processes have ceased in 
the vertical direction. 

(4) Whereas three- and even four-armed modes can dominate at the initial stage of evolution 
of an initially cool disk (Qt{^ < 2L) < 1; see Figs. 2b, 2c), the spiral pattern changes if we 
choose a subcritical initial state (i.e., a state that is unstable, but not so cool, where Qt 
The spatial structure of the perturbations depends to a considerable degree on the relative mass 
of the spheroidal subsystem. Figure 3 shows the distributions of the logarithm of the disk surface 
density at various times. The two-armed mode is dominant in this model, although the m = 3 
harmonic is also important, especially in initial stages of the evolution. It is typical for the spirals 
to join in the outer region of the disk to form a ring-shaped structure. 

(5) Models with low-mass spheroidal subsystems starting from a very cool initial state undergo 
a substantial mass redistribution in the disk, accompanied by a change in the form of the rotation 
curve V(r) in the process of heating and formation of a bar. In this case, the final distribution of 
the surface density a(r) differs strongly from an exponential law (Fig. 4). 5 

Another feature of model disks starting from a very cool initial state (Qt < 1) is that the 
velocity dispersion at the end of the computations (after 10-15 rotations) is somewhat higher 
than is required for gravitational stability. This is due to heating by collective processes — high- 
amplitude wave motions that develop in the presence of strong instability. When the disk is 
heated and reaches marginal stability, these perturbations die out, but the wave-decay process 
has a certain inertia: the velocity dispersion is already high enough to maintain the stability 
of the disk, but the spiral waves have not yet decayed (as is confirmed by Fourier analysis of 
the density perturbations in the disk) and continue to heat the disk. Therefore, to obtain the 
minimum velocity dispersion required for disk stability, we used an iterative algorithm to make 
the initial velocity dispersion approach the stability limit. 

Our iterative approach is based on a series of several (two to four) consecutive computations, 
each of which starts with an initial dispersion that is somewhat closer to the critical level than 
that for the previous computation. For each radius r, we chose an initial velocity dispersion that 
was intermediate between the initial and final values (after five to ten rotations) in the previous 
computation. 

As expected, the stability limit also depends on the initial disk thickness 6 . If the disk is 

5 It is possible, in principle, to choose an initial density profile that yields an exponential profile at the end of 

the computations (dashed line in Fig. 4); however, this approach is rather artificial. It seems that the stellar disks 

of most of the galaxies do not pass through a stage of strong dynamical instability. 

6 The initial thicknesses of the disks of real galaxies depend on the conditions under which they were formed. 

For example, if the collapse of a gaseous disk was accompanied by violent star formation before it reached a quasi- 
stationary state, the resulting stellar disk may be much thicker and have a higher velocity dispersion c z than the 
minimum value required for stability. However, such a scenario is difficult to reconcile with the existence of very 
thin (h z /L < 0.2) stellar disks in edge-on galaxies. The dependence of the disk thickness on the relative mass of 
the halo is also consistent with the assumption that the stellar-velocity dispersion is close to the value expected for 
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initially thick (with the vertical scale height h z >0.2L) and unstable only in its plane (c r < c£ ), 
its heating proceeds more slowly and ceases at lower radial-velocity dispersions than in the case 
of an initially thin disk. This is due to two factors: the stabilizing effect of the finite thickness of 
the disk and the slowness of the relaxation in the vertical direction compared to the heating in 
the plane of the disk. 

The results of dynamical simulations allow to determine the disk parameters at the gravitational- 
stability limit (when the velocity dispersion ceases to change, after 5 4- 20 rotations of the outer 
edge of the disk). 

5.2 Gravitational-Stability Threshold for Bulgeless Models 

We computed a series of bulgeless galaxy models, differing in the relative mass /i = Mh/M^ and 
the halo scale length a, with a fixed radial disk scale L. 

When a >L, the rotation curve has an extended section V c (r) (which we will call the "rigid- 
rotation" section), which makes a transition to a plateau V c ~ const at r >2L (curve 1 in Fig. 5a). 
Figure 5a shows for a bulgeless galaxy with \i = 1 the radial distributions of the circular rotational 
velocity V c (r) (curve 1 ), mean particle rotational velocity V(r) (curve 2), radial-velocity dispersion 
c r (r) (curve 3), and the differential-rotation parameters S = computed separately for V c (r) 

and V(r) (4 and 5). 

Knowing the final velocity dispersions, disk density, and rotational velocity, it is possible to 
compare the Toomre stability parameters Qt (1) for the corresponding model with the analytically 
derived local criteria discussed above and compare these criteria with each other. Note that the 
analytical criteria were derived assuming that the circular rotational velocity V c (r) differs little 
from the mean velocity V(r) of the gravitating point masses. However, the difference between 
these velocities can be quite significant for collisionless disks. We therefore determined the stability 
parameters for both V c (r) and V(r). 

Figure 5b shows for this model the computed Qt and the Toomre parameters required for 
marginal stability according to the criteria of Polyachenko, Polyachenko, and Strel'nikov (PPS) 
Qj 3 ^ (6) and Morozov (5) calculated separately using the circular velocity and the mean 

rotational velocity of the particles. In the case of a massive halo (fi >2), the difference between 
the stability parameters for V c (r) and V(r) is small, but it can become appreciable for /i < 1 
(Figs. 5a, 5b). This is due to the fact that the difference V c — V increases with increasing velocity 
dispersion in accordance with (18). 

The condition S ~ Qj^ <Qt ~Qt'^ is satisfied both in the case considered here and for most 
of the models (recall that Qt = S is the Toomre condition for marginal stability with crude 
allowance for nonaxisymmetric perturbations; see Section 2). 

Figures 5c and 5d compare the stability criteria in a different way. These figures show the 
radial dependences of the ratio of the model velocity dispersion c r of the disk when it becomes 
stable to the critical velocity dispersion for the various criteria considered in Section 2. A perfect 
agreement with the model 

dependences would correspond to the ratio c r /c CT ^ = 1. 

marginal stability [7]. 
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Although none of the criteria explain the model dependences Qrij') at all r, the closest results 
are given by the PPS criterion generalized for the case of finite thickness (2) and the Morozov 
criterion (5) (symbols 15 and 16 in Fig. 5c). When the stability criteria are determined in terms 
of the mean rotational velocity (Fig. 5b), the radial variations of c r /c cr ^ conserve their general 
form (Fig. 5d). Note that for the interval 0.1 < r < 0.8 the ratio c r /c CT ^ is closest to unity for the 
Morozov criterion (curve 16). It is striking that the PPS criterion (curve 12) overestimates c cr j4- in 
all cases, but c r /c cr ^ remains nearly constant over a wide range of r, enabling easy introduction 
of a correction factor for c cr ^. 

The important result is that the radial dependence of the Toomre parameter Qt{t) computed 
for the circular velocity V c (r) shows qualitatively similar behavior for all the bulgeless models we 
considered (Fig. 5e). Moreover, Qrij') remains approximately constant at the level Qt — 1.2 -r 1.6 
in the interval <r/L <2 (see Fig. 5e). At the disk periphery (r >2L), Q T increases with radius, 
reaching Q T ^ 2.5 3 at the edge of the disk (r ~ 4L). However, the scatter of the Qt values for 
the various models is small, enabling us to adopt the following function as a lower envelope 

QP =A + A 1 .(£j+A 2 .{£j 2 , (19) 

where A = 1.25, A 1 = —0.19, and A 2 = 0.134 (bold curve in Fig. 5e). This function reaches its 
minimum, equal to 1.2, at r/L = 0.7. 

The observational estimates of the stellar-velocity dispersions are usually based on data for 
the inner region of the disk r < 2L. Therefore, the approximate constancy of Q T in this region 
makes this quantity a convenient tool for estimating the surface density and mass M d of a disk 
(assuming it is stable) based on the known rotation curve of the galaxy and the radial exponential 
disk scale length L. The equation relating these quantities is 

c r = Q?(r) 3-36GW-r/L) } (2Q) 

where Qj) for galaxies with the extended region of increasing V c (r) is determined by (19). In 
turn, the disk mass estimate M d = IugqL 2 makes it possible to separate out the mass of the dark 
halo within a specified radius. 

5.3 Models with Formation of a Bar 

In the case of a low-mass or very "loose" halo and no bulge, the evolution of an initially cool 
dynamical model inevitably results in the development of a bar. As a rule, the formation of a 
bar involves an appreciable radial redistribution of the mass, whose intensity increases when the 
velocity dispersion in the initial state decreases. The outer boundary of the disk shifts outward. 
As a result, the radial scale for the azimuthally averaged surface density L = — (dlncr/dr)^ 1 varies, 
which can lead to deviations from an exponential profile (when L = const throughout the disk). In 
the innermost region (r <0.3L), the surface density can increase appreciably over its initial value 
(Fig. 4), bringing about a decrease of Qt in accordance with (1). A similar conclusion concerning 
the increase in the central density during the formation of a bar was also reached for models with 
"living" (evolving) haloes [55]. Numerical simulations have shown that, under certain conditions, 
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the interaction between the bar and "living" halo can appreciably affect the dynamical evolution 
of the bar [56] . 

The stability criteria considered above were derived for and applied to axisymmetric disks. 
However, it is of interest to formally compute these criteria based on azimuthally averaged param- 
eters in the absence of axial symmetry in the inner region — that is for barred galaxies. As an 
example, Figure 6 shows the radial distributions of the parameters of a disk that has developed 
a bar (Figs. 6a-6c show results for a model with a low-mass halo and Figs. 6d-6f results for a 
model with a halo whose initial mass within r = 4L is twice that of the disk). The radial distri- 
butions of the surface density in the bar region can differ from those outside the bar (Figs. 6c, 
6f). The development of a bar usually results in additional thickening of the disk. The function 
Qt(j') depends on the bar parameters (essentially, on the initial conditions), however for r >L the 
condition Q T >1.5 is satisfied, and the radial dependence Qt(j') is qualitatively similar to that 
observed for an axisymmetric disk. The development of a bar in models with more massive haloes 
is accompanied by weaker deviations from the initial exponential surface- density profile (Fig. 6f). 

The characteristic dependences for the model with a massive halo (Figs. 6d-6f) differ little 
from the cases shown in Fig. 5 for barless models. As expected, suppressing the bar mode requires 
a higher initial velocity dispersion c r or a more massive halo, in accordance with classical concepts 
(see, e.g., [29, 44]). 

5.4 Models with Bulges 

Let us consider now model galaxies whose rotation curves in the central region (r < L) are deter- 
mined primarily by the bulge. As a rule, the rotations curves outside the bulge (r > L) are almost 
flat (V ~ const). Figure 7 shows typical model radial dependences of the disk parameters for this 
case. In the presence of a bulge, Qt increases strongly in the central region of the disk, where the 
dynamics are determined by the bulge potential (Figs. 7b, 7e). However, outside this region, at 
r ~ (1^-2) x L, the radial dependence of Qt retains its form [see (19), Fig. 5e]. At larger distances 
Qt increases monotonically with r up to values 2.5 -j- 3.5. The resulting radial dependence of Qt 
is typical of systems with bulges that are not too massive (M b /M^ <0.3) and not very extended 
((r 6 )max <L). 

The more massive and more compact the bulge, the higher the value of Qt at the disk center. 
This increase of the Toomre parameter is due mainly to the growth of the epicyclic frequency 
cE. The radial-velocity dispersion in the central region also increases somewhat with increasing 
mass of a strongly concentrated bulge. Additional heating of the center of the disk (r <0.5L) 
in the model shown in Fig. 7a has no direct connection with gravitational instability. At the 
initial stage of its evolution, this model, where fj, — 1, a bar develops, which disappears soon as a 
result of scattering of particles passing near the concentrated nucleus (the latter has a scale length 
b = 0.01). This mechanism for bar disruption is similar to the effect produced by a black hole [57]. 

As a result, the disk undergoes additional heating and becomes thicker in the central re- 
gion (Fig. 7f). The degree of this additional heating and the circular velocity V c at the disk center 
are very sensitive to the bulge parameters, first and foremost, the core radius b (in particular, 
the bar disruption mentioned above does not take place if b >0.05L). Therefore, at r <0.5L the 
Toomre parameter Qt can vary very strongly in different models (Fig. 7f). 
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5.5 Differential Rotation as a Factor Increasing the Threshold for 

Gravitational Stability 

Here we consider two limiting cases: a disk that rotates almost rigidly and a quasi-Keplerian disk 
whose rotational velocity decreases as V oc r~ 1//2 . Since differential rotation is a destabilizing 
factor, rigidly rotating disks are expected to have lower Qt, other conditions being the same. 

Galaxies usually exhibit an extended interval of almost rigid body rotation if the mass of the 
disk dominates throughout most of the disk (r <2L), so that suppression of the bar instability 
requires stronger disk heating than in the presence of massive spherical components. Therefore, to 
elucidate the role of differential rotation, we analyzed a model in which the development of a bar 
is suppressed by a massive halo, yet the outer part of the disk (r > L) contains a rigidly rotating 
section. The differential-rotation parameter in this model S(r > L) = 2VL/de = 1.1 4 1.2 is close 
to S — 1, which corresponds to completely rigid rotation (Fig. 8a). Figure 8a shows the radial 
dependences of Q T , Q^, Qt^ , and S characterizing the stability of the system. It is clear that, 
in general, the disk becomes stable at lower values of Qt (curve 6) than in the cases considered 
above. In this case we have Qt — 1 and Qt — 2 at the disk center and periphery, respectively. 

Consider now another limiting case, corresponding to a nearly Keplerian rotation curve (Fig. 8b). 
Such behavior is very rarely seen in the rotation curves of real galaxies, but we are interested in 
the fundamental question of how strong differential rotation affects the minimum velocity disper- 
sion required for stability. To produce quasi-Keplerian disks, we introduced massive concentrated 
components into our dynamical model. The resulting series of models is a natural extension of 
those with very massive bulges (M b > (2 -4 4)M^). Figure 8b shows the results obtained for one 
of such models. The Toomre parameter increases appreciably for disks with strong differential 
rotation, so that the condition Qt > 2 is necessary for stability, even in the region where Qt is 
minimal (r ~ (1 4- 2)L). 

Since the central region is dominated by the spheroidal system, no bar develops, and there is 
no increase in the velocity dispersion there, which was mentioned in Section 5.4 (Fig. 7a). 

6 DISCUSSION AND CONCLUSIONS 

We have carried out numerical simulations of the evolution of disks that are initially unstable, 
collisionless, and in equilibrium for a whole series of three-dimensional models. The aim of our 
analysis was to compare the minimum radial-velocity dispersions c r at which the disks reach 
their final, quasi-stationary state (as a rule, after 5-10 rotations of the outer edge of the disk) 
with the velocity dispersion Ct = 3.36nGcr/a£ required to suppress gravitational instability of a 
thin disk against axially symmetrical perturbations (the Toomre criterion), as well as with other 
local stability criteria. We separately analyzed the influence of the bulge, dark halo, differential 
rotation, and initial disk parameters on the radial behavior of Qt = c t /ct- 

Our numerical simulations show that differential rotation and inhomogenity of the disk, the 
global nature of perturbations, and the finite thickness of the disk can change the local velocity 
dispersions (or the local Toomre parameters) of marginally stable disks by >50%. The numerical 
models enable to carry out a separate analysis of the effects of each of these factors. 
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The models considered here clearly demonstrate rapid heating of initially unstable disks during 
the formation and destruction of transient spiral arms. The low-contrast "remnants" of these arms 
can be traced over more than ten rotations. A similar character of disk evolution we obtained for 
models which used both different numbers of particles and a different computational algorithm 
(PP method). 

As expected, the minimum radial-velocity dispersion at the end of the simulations (expressed 
in units of the circular velocity) is higher in the models where the relative mass of the halo is lower, 
the initial disk thickness is less, and the degree of differential disk rotation is higher. Although the 
radial dependence of Qt differs for different models and is determined primarily by the relative 
mass and degree of concentration of the spherical components, in all cases, Qt{t) passes through a 
minimum Q T ~ 1.2 4 1.6 at a galactocentric distance of (l-=-2)L, and this behavior depends only 
slightly on the choice of model. This fact can be used to approximately estimate the density (and, 
consequently, the mass) of a galactic disk (or to put limits on these quantities) from observed 
radial- velocity dispersions using equations (19) and (20) without the use of numerical simulations 
or analytical stability criteria. If the halo is not too massive (M^/Mj <2), a disk that begins 
its evolution from a strongly unstable state rather than from a subcritical state can experience 
a substantial mass redistribution over one or two rotations. The fact that galaxies usually have 
exponential brightness distributions indicates that the formation of stellar disks, apparently, was 
not accompanied by the development of strong gravitational instability. 

The thicker is the disk initially, the lower is the minimum radial velocity dispersion c r , which 
determines its stability. Therefore, the minimum critical dispersions in both coordinates z and r 
are reached if the disk begins to evolve from a subcritical state. It's valid both for gravitational 
perturbations in the plane and for bending perturbations (global, primarily axisymmetric mode, 
and small-scale ones). 

The disk evolution traced by the numerical models clearly demonstrates the interrelations of 
processes at different r. Initial stability of the central disk region slows down while strong initial 
instability accelerates the increase of the velocity dispersion at the disk periphery, having, however, 
no significant effect on the final state. 

Existing analytical criteria for the stability of thin disks are based on other considerations, 
and, strictly speaking, there is no reason why they must coincide with the Qt values we have 
found, due to the local nature of the criteria applied, as well as due to the use of two-dimensional 
analytical models, while dynamical models treat the disk heating in three dimensions. The model 
velocity dispersions differ significantly from those predicted by various criteria (Figs. 5b-5d, 6b, 
6e, 7b- 7d, 8a, and 8b). However, it's evident that the formal application of these criteria will yield 
results that are correct to order of magnitude. In particular, like model estimates, they imply that 
Qt > 1 at all r and that this parameter increases at the disk periphery. However, in none of the 
cases considered do the Qt(j') relations derived using local criteria agree with those found in this 
work for all radial distances (0 < r < 4L), and that they can differ from the model estimates by 
more than a factor of 1.5 at some distances r. 

After passing through its minimum, the Toomre parameter Qt{t) for a disk that has reached 
a stable equilibrium state increases monotonically with r up to Qt — 2 4- 3 at a radius of (3 -j- 4)L 
(this radius is usually close to the optical boundary of the disk in real galaxies). Note that a similar 
increase of Qt{t) from 1.5 to 2.5 at the periphery was obtained earlier by Pichon and Lynden- 
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Bell [49] using a different method (by constructing stable equilibrium distribution functions for 
flat systems over a wide range of free parameters) and has also been mentioned for a number of 
dynamical models [5, 21, 24, 44-47]. 

The disks of galaxies whose halo are not too massive (// < 1.5) are subjected to bar instabil- 
ity. The suppression of the development of a long-lived bar in galaxies with low-mass spherical 
components requires not only a higher velocity dispersion (as shown long ago by Ostriker and 
Peebles [29] in their classic work) but also a higher value of the Toomre parameter Qt, compared 
to galaxies possessing a massive halo or bulge. Moreover, the development of a bar is accompa- 
nied by an increase in Qt even in regions of the galaxy beyond the bar. This effect is apparently 
associated with the system being "overheated" above the threshold level: the development of the 
bar is accompanied by a redistribution of mass in the disk, so that the local decrease of the disk 
surface density for the given velocity dispersion results in an additional increase in Q T . However, 
in this case, our models fail to yield bona fide quantitative estimates, since the parameters of 
the bar depend strongly on the initial conditions in the system. We have not analyzed here the 
conditions for the development of the bar mode in detail. 

Our results make it possible to use numerical models to estimate the degree of dynamical 
heating of stellar disks and to derive estimates for the a priori unknown mass and density of a 
galactic disk from the observed velocity dispersion of the old stars making it up. 
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Figurel.Time dependence of the radial-velocity dispersion c r (r) for particles at 20 different 
galactocentric distances (1 is the central zone, and each zone has a width of 0.05 R) at the initial 
stage of heating of a cool, thin disk (t = 10 corresponds to ~ 3 rotations of the outer edge of the 
disk) for two methods of computing the gravitational force with N = 80 000: (a) the PP algorithm 
and (b) the TREE-code algorithm. 
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Figure2. Evolution of an initially unstable thin disk (/j, — 3 ; a — L). (a) Time dependence of 
the radial-velocity dispersion in 20 radial zones, (b) Development of a system of transient spiral 
waves in an initially axisymmetric disk. The positions of point masses in the (x,y) plane at various 
times are presented (since the model consists of 5 x 10 5 particles, only some of the point masses 
are shown), (c) Time dependence of the integrated amplitudes A{m) of the Fourier harmonics for 
various modes m = 2,..., 6. The parameters A reach their maximum values at 1 <t <4 7 and it 
is in this time interval that the perturbations reach their maximum amplitudes (Fig. lb). After 
t > 4, A decreases, which corresponds to the beginning of wave dissipation and the slowing down 
of the growth of c r (Fig. 2a). 
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FigureS. Distribution of the logarithmic surface density log a at various times from t = 1 to 
t = 30 for the initial subcritical state with Qt > 1. 
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Fig 4- Radial dependences of the surface density for an initially cool model disk in which there is 
subsequently a substantial density redistribution. The thin solid curve shows an exponential profile, 
the bold solid curve is the final profile averaged in azimuth after strong heating and development 
of a bar, and the dashed curve is the initial profile that produces a quasi- exponential distribution 
at the end of the computations. 
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Fig. 5. Parameters of the disk that has evolved to a stable state at the stability limit, for bulgeless 

models, (a) Radial dependences of the circular velocity V c (r) (curve 1), disk rotational velocity 

V(r) (curve 2), radial-velocity dispersion c r (r) (curve 3), and parameters S = 2Q/<£ computed 

for the circular velocity ( curve 4 ) and disk rotational velocity ( curve 5 ), respectively. Results are 

shown for the model in which the mass of the halo within r < 1 = 4L is equal to that of the disk 

and the halo scale length a is equal to the radial disk scale length L. (b) Radial dependences of 

the Toomre parameter at the stability limit computed using various stability criteria for the model 

(p) 

shown in Fig. 5a: 6 Qt determined from V c , 7 Qt determined from V , 8 Q T determined from 
V c , 9 Qj?^ determined from V, 10 Qj^* determined from V c , and 11 Q^* determined from V . 
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Fig.5(cont.)(c) Radial distributions of the ratio c r /c cr ^ of the radial-velocity dispersions in the 
numerical models (curve 3 in Fig. 5a) to the critical values computed using the circular rotational 
velocity V c (r) (curve 1 in Fig. 5a) for 12 the Polyachenko-Polyachenko-Strel'nikov criterion (6) 
and 13 the Morozov criterion (5). Also shown are 14 c r / \ct1VL / 'as) derived using the simplified 
criterion (4); 15 c r normalized to c^Q^, which is a generalization of criterion (6) for the case 
of a finite-thickness disk, derived using (2); and 16 the ratio obtained for the Morozov criterion 
(5) in the case of an infinitely thin disk (D = 1). (d) Same as Fig. 5c using the disk rotational 
velocity V(r) (see curve 2 in Fig. 5a). 
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Fig.5(cont.) (e) The Toomre parameter Qt{ t ) a t the stability limit for a series of bulgeless 
models with various halo parameters. The bold thick curve is computed using (19). 
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Fig. 6. Bulgeless models that develop bars in the course of their evolution, (a)-(c) = 0.7Mj 7 
and a = 2L; (d)-(f) = 2Mj and a = 1.6L. Same notation as in Fig. 5. Plots c and f show 
azimuthally averaged surface- density profiles at the initial time (dashed) and after the development 
of a quasi-stationary state with a bar at t > 25 (solid). 
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Figurel. Disk parameters at the stability limit for models with bulges, (a)-(e) Radial de- 
pendences of the disk parameters (same notation as in Fig. 5). The mass of the halo within 
r < 1 = 4L is equal to that of the disk, and the halo scale length is a = 3.6L. The bulge param- 
eters are Mb = 0.24M^ 7 b = 0.04L, and {n)max = 0.5L. (f) Toomre parameter Qr(j') computed 
using V c (r) at the stability limit for the series of models with bulges. The bold solid curve is based 
on (19). (g) Edge-on view of the disk at the end of the computation of the model shown in a. The 
dots show the positions of the particles. The model develops an appreciable bulge-like feature in 
the central region of the disk. 
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Fig. 8 Radial dependences of the disk parameters in the adopted units for the cases when (a) 
a significant part of the disk rotates in accordance with a quasi-rigid law, and (b ) the rotation at 
r > L = 0.25 is close to Keplerian: 1 the circular velocity V c (r); 2 the rotational velocity of the 
disk (stars) V(r); 3 the radial-velocity dispersion c r (r); 4 the parameter S computed using (3) 
for the circular velocity; 6 the Toomre parameter Qt computed using the circular velocity; 8 the 
critical Toomre parameter computed using (6) with V c (r); 10 the critical Toomre parameter 
computed using (5) with V c (r). The notation is analogous to that in Fig. 5. 



